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QUASI-LINEAR PARAMETER VARYING 
REPRESENTATION OF GENERAL AIRCRAFT DYNAMICS 
OVER A NON-TRIM REGION * 

Jong-Yeob Shin 1 


ABSTRACT 

To apply linear parameter varying (LPV) control synthesis and analysis to a nonlinear 
system, it is required that a nonlinear system be represented in the form of an LPV 
model. In this paper, a new representation method is developed for the construction 
of an LPV model from a nonlinear mathematical model without the restriction that 
an operating point must be in the neighborhood of equilibrium points. An LPV model 
constructed by the new method preserves local stabilities of the original nonlinear 
system at “frozen” scheduling parameters and also represents the original nonlinear 
dynamics over a non-trim region. An LPV model of the motion of FASER (Free-flying 
Aircraft for Subscale Experimental Research) is constructed by the new method. 

NOMENCLATURE 

physical parameters 

u, v, w : X-, Y-, Z-components of aircraft velocity in body axes 

p, q, r : X-, Y-, Z-components of aircraft angular velocity in body axes 

4>, 9, if} : Position angles 

a, (3 : Angle of attack, side-slip angle 

m, J, g : Total mass, moment of inertia matrix, gravity constant 

V, T : Total speed, thrust force 

L, M, N : X-, Y-, Z-component of moment 

control surfaces 

5 e : elevator deflection angle 

5 r : rudder deflection angle 

5 ra , 8 ia '■ right and left side aileron deflection angle 

S r f, 5if : right and left side flap deflection angle 

This work was supported by the National Aeronautics and Space Administration under NASA Cooper- 
ative Agreement NCC1-02043. 
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aerodynamic coefficients 

Cx, Cy, Cz : X-, Y-, and Z-component force aerodynamic coefficients in body axis 

Ci, C m , C n : Roll, pitch, and yaw moment aerodynamic coefficients 

Cx 0 , Cz 0 , C mo : Aerodynamic coefficients for X- and Z-forces and pitch moment 

Cx q , Cx Se , Cx Sra '■ X-force aerodynamic derivatives due to q, 5 e , and 5 ra 

Cx s j Cx s , Cx s : X-force aerodynamic derivatives due to Si a , 5 r f, and 5[f 

Cy 0 , Cy p , Cy r : Y-force aerodynamic derivatives due to /3, p, and r 

Cy Sr , Cy Sra , Cy Si : Y-force aerodynamic derivatives due to 5 r , 5 ra , and 5i a 

Cy s , Cy s : Y-force aerodynamic derivatives due to 5 r f and 5if 

Cz q , Cz Se , Cz Sra : Z-force aerodynamic derivatives due to q, S e , and 5 ra 

Cz s j Cz s , Cz s : Z-force aerodynamic derivatives due to 5i a , 5 r f, and <5// 

Ci 0 , Ci p , Ci r : Roll moment aerodynamic derivatives due to (3, p, and r 

c kr , Ci Sra , C h : Roll moment aerodynamic derivatives due to 5 r , S ra , and 5{ a 

Ci s , Ci s : Roll moment aerodynamic derivatives due to 5 r f and 5if 

C mq , C mge , C mSra '■ Pitch moment aerodynamic derivatives due to q, S e , and 5 ra 

Cm 5l > C ms , C mg : Pitch moment aerodynamic derivatives due to 5i a , 5 r f, and 5if 

C ri0 , C np , C nr : Yaw moment aerodynamic derivatives due to (3 , p, and r 

C nSr , C ng , C ndla : Yaw moment aerodynamic derivatives due to S r , S ra , and 5( a 

C ns , C nSlf ■ Yaw moment aerodynamic derivatives due to S r f and 5if 

1 INTRODUCTION 

Gain-scheduled control techniques have been applied for control of a nonlinear system over 
decades [22, 28, 12]. One of the gain-scheduled control techniques promising global stability 
of a controlled system is an LPV control synthesis methodology [2, 6, 32, 25]. An LPV 
control methodology is also particularly appealing in that robust linear control techniques 
in an LPV system can be applied to a nonlinear system using well-developed linear robust 
control concepts [34] without the loss of stability guarantees. To apply an LPV synthesis and 
analysis methodology to a nonlinear system, it is required that its nonlinearity be represented 
in the LPV form over the entire operating envelope. 

A few studies [2, 30, 24, 13, 15] have been reported considering LPV representations of 
nonlinear systems. In Refs. [2, 14, 21], an affine-LPV model has been introduced, whose 
system matrices are affine functions of parameters. To represent nonlinear dynamics as an 
affine-LPV model, the convex set of a constructed LPV model is generally larger than actual 
nonlinear dynamics. This may lead to conservative results in control synthesis. An alterna- 
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tive method of LPV representation is a grid-LPV model whose system matrices are functions 
of the parameters at all grid points over the entire parameter space [6, 32, 25]. This LPV rep- 
resentation has been used for a grid-LPV synthesis method in which linear matrix inequality 
constraints are evaluated at all grid points. Generally, there is a trade-off between accuracy 
in representing the original nonlinear dynamics and computational cost (memory and com- 
putational time) for a grid LPV representation method. The grid-LPV representation and 
control synthesis methods have been used for several control problems: pitch-axis missile 
autopilots [33, 24], F-14 aircraft lateral-directional axis during powered approach [5, 3], and 
turbofan engines [31, 4], In this paper, we will focus on a grid-LPV representation of a non- 
linear system since the representation can be applied to a large class of nonlinear systems 
without the restriction of affine-function dependence of scheduling parameters. 

There are a few methods [10, 24, 30] to generate a grid-LPV model from a nonlinear 
mathematical model. The most common way to generate a grid-LPV model is the first 
order approximation of a nonlinear system around equilibrium points. This is called the 
“Jacobian linearization method” which has been used in many applications [10, 20, 25, 15]. 
As noted in Refs. [10, 20, 25, 15], an LPV model generated by the Jacobian linearization 
method can capture nonlinear dynamics only around equilibrium points. Another method 
is a “state transformation” method which has been introduced in Ref. [24], In this method, 
state coordinates are changed along with an equilibrium trajectory of scheduling parameters 
to remove nonlinear dependence on states and inputs of a nonlinear system. The method 
can generate an accurate LPV model (not using the first order approximation) of a nonlinear 
system. The disadvantage of the method, however, is that it cannot generate an LPV model 
outside the equilibrium manifold. 

One approach to develop an LPV model over the entire operating envelope including 
outside the equilibrium manifold has been introduced in Refs. [30, 26, 15]. It is called a 
“function substitution” method in which a nonlinear mathematical model is replaced by a 
linear combination of functions of scheduling parameters at all grid points over the entire 
operating envelope. The linear-combination functions are then reformulated into LPV form. 
To decompose nonlinear functions, a single reference point (one of the equilibrium points) 
is used, as shown in Refs. [26, 15]. The LPV model generated by the method in Refs. [26, 15] 
is strongly dependent on the single reference point. Also, the constructed LPV model may 
not capture the local stability of the original nonlinear system at an equilibrium point. In 
this paper, we develop a new function substitution method to remove single reference point 
dependence on an LPV model and to preserve the local stability of the original nonlinear 
dynamics, while representing the nonlinear dynamics over the entire operating envelope. To 
preserve local stability, an LPV representation is formulated as a linear optimization problem 
which is solved by using the linear matrix inequality (LMI) solver interface “YALMIP” [11], 

The remainder of the paper is organized as follows: in Section 2, the class of nonlinear 
systems and a quasi-LPV model are defined and the preliminary background is explained. 
In Section 3, the current LPV representation methodologies are described in order to math- 
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ematically define the problem of interest. In Section 4, a new method is described to resolve 
the problem. In Section 5, the nonlinear dynamics of FASER is described and the suggested 
method is applied for constructing an LPV model of the longitudinal/lateral-directional 
motions of the aircraft. The paper concludes with a brief summary in Section 6. 

2 PRELIMINARY 

General Nonlinear Aircraft System 

In this paper, general aircraft dynamics are represented by the following nonlinear mathe- 
matical model. 

x — F(xi)x + G(xi)u + h(xi) (1) 

with a state vector x G 7 Z Hx , and a control vector u G 7 Z Uu . F(x i), G(x i), and h(x i) 
are continuous mapping functions: TZ nx i i— > 7£ n:tXnx , JZ n ^i i— > pi nxXnx , and 1Z nx i > VJ ' X , 
respectively. In Eq. (1), a state vector x is written as x T = [xf xlf]. Note that the 
subvector x\ is part of the state vector and also a scheduling parameter vector. Assume that 
the scheduling parameter variations are in the compact set V in FJ lx \ , which is called the 
operating envelope, hereafter. In the operating envelope, the equilibrium locus is defined as 
the set of equilibrium points (aq, x 2 , u) such that 

0 = F(xi) ^ 1 + G(xi)u + h{£i). (2) 

pi_ 

Definition 1 (local stability): The linearized system of a nonlinear system in Eq.(l) is 

x = Ax + Bu, (3) 

where matrices A and B are obtained by Jacobian linearization about an equilibrium point. 
When the real parts of all eigenvalues of A are negative, the nonlinear system is locally stable 
about the point. 

Even though this is a very common definition, it is described here to emphasize the rela- 
tion between time responses of a nonlinear system with small perturbations on states and 
the stability results based on the sign of all eigenvalues. Suppose the linearized system at 
an equilibrium point indicates stability by definition 1. It is obvious that nonlinear time 
responses of the system with small perturbations on states with input u = u (equilibrium 
value) should converge to the equilibrium point. In practice, it is possible to test local sta- 
bility based on time responses unless the system dynamics are very slow. Note that the 
local stability applies for the neighborhood of the equilibrium point only. To determine local 
stability for arbitrary points over the entire envelope, time responses of a nonlinear system 
are simulated with perturbations on the state. This is referred to in this paper as “extended 
local stability” , and is defined as follows. 
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Definition 2 (extended local stability) '.Suppose there exists the set £ of pre- calculated 
equilibrium points. Given any point X\ in the operating envelope V , there exists an equilib- 
rium point X\ r such that 


X\ r — argmin £ie£ \\W[xi - Xi]||, Xi G V. (4) 

where W is a weighting matrix to compensate physical unit difference of each state. The 
trajectory of the nonlinear system in Eq. (1) is calculated with the initial condition ofx( O) 2 = 
[xf xff r } mid the control input u = u r . Here, equilibrium values of x 2r and u r are associated 
with xi r . Suppose there exists small number e such that 


x(t) 



t > T > 0, 


( 5 ) 


then, the system is stable with respect to the point (xi r , x 2r ) around the neighborhood of the 
point xi in the operating envelope. 


Note that the extended local stability about a grid point does not represent global stability 
of the nonlinear system [9] and is also strongly dependent on both x\ T and X\. When the 
time history of x(t) converges to another point instead of {x\ r , x 2r ), stability is not captured 
by definition 2. From the extended local stability definition, we can only determine whether 
the system is stable around the grid point. It is noticed that the minimum solution x\ r may 
not be unique with respect to X\. This implies that the result of the extended local stability 
may not be unique. 

2.1 Quasi-LPV System 

In this subsection, a quasi-LPV system is briefly described to carry out the LPV repre- 
sentation of the class of nonlinear systems in Eq. (1). A quasi-LPV system is defined as 
follows: 

Definition 3 (Quasi-LPV system): A quasi-LPV system is defined as 


X 


A(xi) 

B( Xl ) 

X 

v . 


C(x i) 

D(x i). 

u_ 


with states x 1 = [xf xf] . Note that the states X\ are scheduling parameters whose variations 
are in the compact set V G 1Z nx i . 

The detailed description of a quasi-LPV system is available in Refs. [24, 20]. When a quasi- 
LPV model is used for analysis or control synthesis for a nonlinear system, it is assumed that 
scheduling parameters can independently vary even though the schedule parameter vector is 
the part of the state [24, 20]. 
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Stability analysis on a nonlinear system can be formulated into a linear matrix inequality 
(LMI) optimization using a quasi-LPV representation of the nonlinear model. When a quasi- 
LPV model in Eq. (6) represents the original nonlinear dynamics, the stability analysis 
problem is described as follows: 

Lemma 1 (global stability) [6] Suppose there exists a positive definite matrix P such 
that 

a t (x 1 )p + pa(x 1 ) < 0 , y Xl ev. (7) 

Then, the system is globally stable over the operating envelope V . 

This is easily proved using Lyapunov function V = x T Px. In practice, the LMI constraints 
in Eq. (7) are evaluated in all grid points over the envelope V. Thus, the frozen dynamics 
at a fixed parameter plays an important role with determining local and global stability of 
the original nonlinear dynamics. 


3 PROBLEM STATEMENT 


3.1 Set of Linearized Models 


Recall that a set of linearized models of a nonlinear system is considered as a quasi-LPV 
model for the control synthesis [25]. One approach to construct a linearized model set is to 
linearize a nonlinear system about equilibrium points using a Jacobian linearization method. 
The set M.jl of linearized models of a nonlinear system Eq.(l) is 


Mjl = {x = 


,d(F(xi)x) 

dx 


, { dU{xi)u i 

dh[x\) 

\x 1 0 Lr, , u - I 

OX 

dx 


-~)x+G(xi)u : x = x—x,u = u—ii} (8) 


where ( x,u ) is an equilibrium point. When scheduling parameters X\ are chosen along the 
locus of equilibrium points, X 2 and u are associated with X\. Then the set of linearized 
models can be written in a quasi-LPV form: 


x = A{x\)x + B{x\)u 


( 9 ) 


where 


A(xi) 


d(F(x\)x) I | dG(x\)u\ , dh(x!) i 

dx \x ' Q x \x,u \ g x \xi 


B(x 1 ) = G'(xi). 


Note that this quasi-LPV model is a collection of linearized dynamic systems around the 
equilibrium points/trajectory [13] since the state and input definitions are changed according 
to an equilibrium point associated with the current scheduling parameters (£i). Advantages 
of this well-known method are: 1) applicable to a wide-class of nonlinear systems without 
restriction of affine dependence of input u on a nonlinear system, and 2) able to represent 
local stability of the nonlinear system preserving the stability of a quasi-LPV system at 
“frozen” parameters described in Lemma 1. 
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The quasi-LPV model constructed by this method can represent the true nonlinear dy- 
namics at the neighborhood of equilibrium points only. Choosing the number of grid points 
along the locus of equilibrium points also plays an important role in the representation of 
the nonlinear dynamics. The optimal way to choose grid points is still unknown. 

3.2 Function Substitution Quasi-LPV Model 

The function substitution method has been introduced in Refs. [30, 25, 15] to construct a 
quasi-LPV model over the entire operating region, including a non-trim region. To perform 
the substitution method, choose an equilibrium point as a reference point (x lr , x 2r , u r ) and 
transform the states into 


Xi — Xi — Xi r , X 2 = X 2 — X 2r , U = U — U r . 

Using Eq. (10), the nonlinear system of Eq. (1) can be rewritten as 


X\ 

x 2 


where 


f( x i) = 


A u{xi) Ai 2 (xi) 

X \ 

+ 

Bi(xi) 

_ A 2 l(xi) A 22 (x\)_ 

X 2 _ 


B 2 (x i)_ 

^11(^1) Ai 2 (x\) 

X\ r 

+ 

Bi(xi) 

^21(^1) ^22(^1). 

_% 2 r _ 


B 2 (x i)_ 


u + f(x i) 


u r + h(x i). 


( 10 ) 


(ii) 


( 12 ) 


The main goal of this method is to reformulate the function f(x\) in Eq.(12) into quasi-LPV 
functional form such that 


f(x i) = E(xi)xi 


ei.(xi) 


Xi 


e n -( x i) 


(13) 


where the matrix E is an unknown matrix to be determined and e t . is the i-th row vector 
of the matrix E. There are infinite numbers of possible solutions to Eq. (13) since this is 
an under-determined problem. In Refs. [25, 15], the matrix E is calculated to minimize the 
variations of each matrix element over the entire operating envelope. 


The final quasi-LPV form of the nonlinear system is 


where 


Xi 

x 2 . 


A f 


Xi 

x 2. 


+ 


B\{x\) 

B 2 M 


u. 


A f = 


A ll(xi) 

.^21(^1) 


2112(2^1) 

^22(^1) 


+ [E I 0 nxn:l;2 _ • 


(14) 


(15) 
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The solution of the quasi-LPV system is close to the solution of the nonlinear system since 
the equality constraint in Eq. (13) is satisfied at all grid points. This is a strong advantage 
of this method. 

This method, however, has disadvantages such as 1) a strong dependence of the quasi- 
LPV model on the reference point and 2) possible misrepresentation of local stability of the 
original nonlinear system over the locus of equilibrium points. The matrix Af of the quasi- 
LPV model can change when a different equilibrium point is chosen as a reference point. 
The stability analysis described in Lemma 1 may fail even though the original nonlinear 
dynamics are exponentially stable. Since all LMI constraints in Lemma 1 are evaluated 
at all grid points, the matrix Af at each grid point must be stable to represent the stable 
nonlinear dynamics. The current limitations of the method in Refs. [26, 15] are choosing the 
correct reference point and representing the local stability of the original nonlinear system 
at a grid point. In the next section, the limitations will be relaxed by a new function 
substitution method. 

4 QUASI-LPV REPRESENTATION DEVELOPMENT 

A new method is described in this paper to construct a quasi-LPV model of a nonlinear 
system over the entire operating envelope including a non-trim region. The model can 
also capture the local stability of the original nonlinear system. In the new method, many 
reference points are chosen along the locus of equilibrium points to preserve the local stability 
of the system. The details are described in the four step process below: 

First, grid points are generated over the entire operating envelope V and each grid point 
is assigned with its own reference point. This allows us to reformulate the nonlinear function 
f(x i) in Eq.(12) into quasi-LPV form. A reference point for each grid point is determined 
by Eq. (4). 

Second, E in Eq. (13) is written in terms of the variables 0* to satisfy the equality con- 
straint. When a grid point x\ is not equal to a reference point aq r , Eq. (13) is rewritten 
as: 



The possible solutions of ej. in Eq. (16) can be written as 

e l = e l P + •A/'Qu - Ur] T )0i (17) 

where the particular solution ef. p is 

e l v = [U ~ %lr\(,[Xl ~ X lr } T [Xi - Xir^fj,, 

8 


(18) 



Af([xi — x\ r ] T ) is the null space of [aq — xi r ] T , and 0,; is a unknown vector associated with 
the null space dimension. Using Eq. (17), the matrix Af is rewritten as: 

A f = A p + [d> T A7 T ([xi - x lr ] T ) j 0 BXn J (19) 


where 


A p — 


An 

A 21 


A12 

A22. 


+ [E p \0 r 


( 20 ) 


the matrix E p is the collection of ef. p and the matrix $ is the collection of 0j. 


When a grid point is equal to a reference point X\ r , notice that f[x 1 ) and x\ are zero. The 
particular solution of Eq. (18) cannot be defined since it is a singular point. For this case, 
the matrix E dehned in Eq. (13) is defined as 


E 


lim 


f(x 1 ) 


xi—>xi r X\ — Xi r 


( 21 ) 


Using Eqs. (15) and (21), the matrix Af is calculated for a quasi-LPV model at the grid 
point. Noticed that the calculated Af is exactly the same as the matrix A obtained by 
Jacobian linearization at the grid point. This is proved after algebraic manipulations with 
Eqs. (11) and (12). 


Third, a linear matrix inequality is constructed to represent stability of the original non- 
linear dynamics with prior information of the extended local stability of the system at every 
grid point. When the original nonlinear system is stable, the following LMI must be satisfied 
for all grid points. 

(A p + [<& T Af 1 | 0 nxn JfP + P(A p + [<h T A f T | 0 nX n X2 }) <0, P > 0, aq G V (22) 

where P is a constant positive matrix. Note that Eq. (22) has unknowns P and in 
multiple form. This makes the LMI constraints not convex in P and <f>. ft is a bilinear 
matrix inequality (BMf) problem. Generally, it is hard to find a global solution to a BM1 
problem. In this paper, the feasible solution of <f> is calculated using the Ruth-Hurwitz 
condition. Then, the LMI optimization Eq. (22) is solvable with the calculated feasible 
solution <f>. The detailed method is described in the next section. 


Fourth, for smoothness of the matrix Af elements, a linear optimization problem is formu- 
lated to penalize the derivatives of the matrix elements due to X\ . The details of constructing 
a linear optimization problem is described in the next section. 


5 EXAMPLE 


In this section, it is demonstrated how to construct a quasi-LPV model using the proposed 
method. The nonlinear dynamics of FASER are represented by a quasi-LPV model. FASER 
(see Fig. 1), developed at NASA Langley Research Center, has conventional high wing and 
tail configuration with 7 ft wingspan [17]. The detailed shape and physical quantities are 
available in Refs. [16, 17]. 
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5.1 Nonlinear Equations of Motion 


In this subsection, general nonlinear equations of aircraft motion [7] are reformulated into 
matrix form described in Eq. (1). 


Force equations of an aircraft are 


u = rV sin (3 — qV sin a cos (3 — g sin 9 + 


qSC x + T 


m 


v = pV sin a cos / 3 — rV cos a cos (3 + g cos 6 sin 0 + 


qSC Y 


m 


w = qV cos a cos (3 — pV sin (3 + g cos 6 cos 0 + 


qSC z 


m 


Kinematic equations of position angles are 


and moment equations are 


p 

q 

= J' 1 

L 

M 

+ J” 1 

r 


N 



0 = p + (q sin 0 + r cos 0) tan 6 

6 = q cos 0 — r sin 0 

ip = (q sin 0 + r cos 0) / cos 6 


JzxQ 3 X yT JyzQ Jzz T JyyQ 3y Z V 


jiryV Jyi/1 J xxV JxyQ JyzP JzxQ 

where the inverse of the inertial moment 3 is 



P 


q 


r 


3~ l = ^~ 


3. 


D 


J J - J 2 

Jyy^zz J y Z 
JxyJzz 3 X z Jyz 


JxzJyz T JzzJxy JxyJyz T 3 xz Jyy 
J J - J ' 2 

XX u 7.7. w 


JxyJyz~ I - JxzJyy JxxJyz JxyJxz 


JxxJyz “ 1 “ JxyJxz 
J J - J 2 

u xx u yy u X y . 


and the scalar Jd is 

Jd JxxJyyJzz 


J J 2 _ J J 2 

*3 xx ^ yz u zz*J X y 


2T T J — T J 2 

^^xy^xz^yz J yy u xz • 


(23) 

(24) 

(25) 


(26) 

(27) 

(28) 


(29) 


(30) 


(31) 


The roll, pitch and yaw moments (L, M, N) are functions of dynamic pressure q, reference 
area S , wing span b, mean chord length c, and moment aerodynamic coefficients Q, C m and 


C n . 


L = qSbCi , M = qScC m , N = qSbC n 


(32) 


After simple algebraic manipulation of Eqs. (23)-(25), the force equations in Eqs. (23)-(25) 
can be rewritten in terms of total velocity V, angle of attack a, and side-slip angle (3 . This 
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leads to 


V 


= — {(qSC x + T) cos a cos / 3 + qSCy sin (3 + qSCz sin a cos (3} 
m 

+ g { cos 9 sin 0 sin f3 — sin 9 cos a cos / 3 + cos 9 cos 0 sin a cos (3) 


( 33 ) 


cos a sin/3 sin a sin ft qSCz cos a qSC x + T 

a = q — p ^ — r H — — sin a 
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V cos (3 


cos / 3 cos / 3 mV cos (3 

cos a cos 9 cos (f) + sin 9 sin a) 


mV t cos (3 


(34) 


h . , qSC Y a qSC x + T . a T,SC Z . a . 

/ 3 — p sin a — r cos a H — cos (3 — cos a sin p — sin p sin a 


mV* 


mV* 


mV 
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— (cos 9 sin (f) cos (3 + sin / 3 cos 9 cos a — cos 9 cos 0 sin a sin (3). 
Vt 


(35) 


From Eqs.(33)-(35), (26)-(29) notice that the states are V, a, (3, body angular rates p, q, r, 
and position angles 0, 9, and 0 for the nonlinear dynamics. 


5.2 Aerodynamic Coefficients 

Force and moment equations in Eqs. ( 23 )-( 25 ), ( 29 ) are dependent on aerodynamic coef- 
ficients Cxi Cy , Cz, Ci, C m , and C n . In this example, the aerodynamic coefficients are 
functions of angle of attack, and control surface deflections such as elevator deflection 5 e 
(deg), right aileron deflection 5 ra (deg), left aileron deflection 5i a (deg), right flap deflection 
S r f (deg), left flap deflection 5if (deg), and rudder deflection 5 r (deg). The aerodynamic 
coefficients are 

C x =C Xo {oi) + C Xq (a)—q + Cx 5e (a) A + Cx Sra (o.)8 ra + C x&la [p)8\ a 

+ Cx Srf ( a)8 rf + Cx hf (p)8i f ( 36 ) 

C Y =Cy 0 (a)(3 + —(C Yp (a)p + C Yr (a)r ) + Cy Sr (a)S r + Cy Sra (a)5 ro + C Ysia {oi)8i a 

+ C Ys rf ( a )8rf + Cy hf ( a)5if ( 37 ) 

C z =C Zo (a) + C Zq (a)—q + Cz Se (a) A + Cz Sra (a)8 ra + Cz ha {a)5 ia 

+ Cz Srf (a)5 rf + Cz hf (a)5 if ( 38 ) 
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( 39 ) 


Ci — Ci p {a)(3 + — (Ci p (a)p + C ir (a)r ) + Ci Sr (a)S r + Ci Sra (a)S ra + Ci ha (a)5i a 
+ Ci Srf {a)8 r f + C iSif ( a)5 t f 

C m =C mo (a ) + (a) —q + C mSe (a)8 e + C mi5ro (Q;)(5 m + C ms ^(a)5i a 

+ C m5 (a)^r/ + C ms ( a)5if (40) 

c » + + CMr) + + + c ^ a)Sl ° 

+ Cn 5rf (<x)$rf + Cn hf ( Ot)5if (41) 

The aerodynamic coefficients Cx a , Cz 0 , and C mo and aerodynamic derivatives are func- 
tions of angle of attack and were calculated by Paul Chwalowski[19] using CMARC [1] and 
VORVIEW [8, 27]. The CMARC is a commercial version of potential flow method developed 
at NASA Ames. The calculated aerodynamic coefficients are compared with the results of 
wind tunnel experiments[17] done at the NASA Langley 12-foot Low Speed Wind Tunnel. In 
this paper, the calculated aerodynamic coefficients are used to develop a quasi-LPV model 
of the aircraft since the comparison results [19] show that the CFD calculation accurately 
predicts the aerodynamic coefficients of the FASER aircraft. 

5.3 Reformulated Nonlinear Equations of Motion 

In this subsection, the nonlinear equations described in the previous section are written in 
matrix form of Eq. (1) to formulate quasi-LPV representation. 

When states X st — [V a (3 p q r 9 0] T and controls U — [T 8 e 5 ra 5i a 8 r f $if 8 r ] T , 
the nonlinear equations of motion can be rewritten as 

X H = A(V, a , /3, p, q, r, 9, 0)X st + B(V, a,(3)U+ h(V, a, (5) (42) 

where the detailed elements of the matrices A, B, and h are written in Appendix A. 

Note that it is possible to generate a quasi-LPV model for the full 6 DOF motion but 
scheduling parameters in a quasi-LPV model would be all states. When all states are schedul- 
ing parameters, a quasi-LPV model can lead very conservative results in control synthesis 
since scheduling parameters are considered independent on each other in LPV control syn- 
thesis [29, 23]. 

To generate a practically useful quasi-LPV model for control synthesis, the nonlinear 
dynamics are decoupled into the longitudinal motion and the lateral- directional motion. 
For the longitudinal motion of the aircraft, state Won is set as Won = [V a q 6] 1 
and control Ui on is set as Ui on = [T 5 e S as 8 fa ] T , where S as = 0.5(<5 m + 8 [a ) and S fs = 
0.5 (S r f + 5if). Positive deflection of all control surfaces except the rudder is defined as trailing 
edge down. For control of the longitudinal motion, the right/left ailerons are assumed to 
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move symmetrically and so do the right/left flaps. Since the other variables are set as zero 
for the longitudinal motion, the nonlinear equation of longitudinal motion is 

X-lon Ai on (V, a)Xi on T Bi on (y, (XjUion 4- h/ on ( V*, Ot) , (43) 

where the elements of the matrices Ai on , Bi on , and hi on are given in Appendix B. 

For lateral- directional motion, the states are Xi at = [/3 p r (j)] T and the controls are 
Uiat = [5a as S fas 5 r ] T , where 5 aas = 0.5(4* - 8 ra ) and 5 fas = 0.5 (5 if -8 rf ). The right and left 
ailerons are assumed to move asymmetrically for control of lateral- directional motion. The 
nonlinear equations of lateral-directional motion is 

Xiat — Ai at (V, a, (5, 4>)Xi at + Bi at {V , a)Ui at (44) 

where the elements of the matrices Ai at and Bi at are given in Appendix B. From Eq. (44), 
notice that the matrix form of the lateral-directional motion is already in quasi-LPV form 
with exogenous scheduling parameters V and a and state-scheduling parameter f3 and 

5.4 Quasi-LPV Model of Longitudinal Motion 

To construct a quasi-LPV model of the longitudinal motion by using the new method, the 
locus of equilibrium points is calculated over the flight envelope defined as 

50 < V < 200 ft/s, 0.5 < a < 10 deg. (45) 

The chosen equilibrium points are numbered as shown in Fig. 2 while the grid points are 
defined as shown in Fig. 2. Using Eq. (4), every grid point is assigned to an equilibrium 
point. The assigned equilibrium point is called a reference point for the grid point. 

To apply the method, the stability information of the nonlinear longitudinal dynamics 
is required. The local stability of the motion around equilibrium point is easily tested by 
checking eigenvalues of the matrix A obtained by the Jacobian linearization method. The 
extended local stability at every grid point can be tested by using nonlinear simulations at 
the initial values of {x\,X 2 ) using the constant input associated with the reference point. 
Based on the stability test results, we conclude that the longitudinal dynamics are stable for 
all the grid points. 

Now, the decomposed function f(xf) of Eq. (12) is evaluated at the grid point X\ using 
Eq. (17), the particular solutions are calculated at all grid points except the equilibrium 
points. Note that the matrix Af of the quasi-LPV model at the numbered equilibrium 
points is calculated by using the Jacobian linearization method. It is required to find feasible 
solutions of <h T at a grid point to solve the BMI problem. We first try to find the feasible 
solution by assigning 4> r as a zero matrix. Note that the particular solution of A p in Eq. (20) 
can be the system matrix of a quasi-LPV model. The eigenvalues of the matrix A p at every 
grid point are shown in the left plot of Fig. 3. It is observed from Fig. 3 that there are 
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unstable eigenvalues of A p at some grid points. Thus, this solution is not appropriate to 
represent the stable nonlinear dynamics of FASER. 

Since the matrix Af in Eq. (19) can be explicitly written as <f> r , the Ruth-Hurwitz stabil- 
ity [18] of the 4th-order polynomial can be used to find feasible solutions of <f> r . The stability 
conditions of the 4th-order polynomial s 4 + as 3 + bs 2 + cs + d are a > 0, b > 0, c > 0, d > 
0, ab > c and dbc — c 2 > da 2 . To use the stability condition, the matrix Af is rewritten as: 


Af(V.a) 


(V,a) 

+ (pi (V, a)ni(y, a) 

a p 12 (V,a) 

+ (pi (V, a)n 2 (V, a) 

a p 13 (V,a) 

a p u (V, a) 

(V,a) 

+ 02 (V, a)ni(V, a ) 

<4 (V,a) 

+ <p 2 (V, a)n 2 (V, a) 

a 23 (V,a ) 

«24 (Via) 

(V,a) 

+ 

CO 

77 

s 

77 

CO’tS 

to 

77 

p 

+ 03 (V, a)n 2 (V, a) 

a| 3 (R,a) 

0 


0 


0 

1 

0 






(46)' 


where 0 i, 02 , and 0 3 are unknown parameters to be determined, ri\ and ri 2 are the elements 
of the null space and ab- is the element of the matrix A p . For convenience, the dependence 
of V and a on parameters is omitted hereafter. 


After some algebraic manipulations, the conditions a > 0 , 6 > 0 , c > 0 , d > 0 can be 
rewritten as 

9a x > °> 9b x > 9l x > 0 , g d x > 0 ( 47 ) 

where x 1 = [1 (pi 0 2 03 ] and the coefficient vectors gj, g gj, and gj are written in terms 
of ab- , ni and ri 2 ■ Here, the detailed expression of the coefficients vectors are in Appendix 
C. The condition ab > c is relaxed by the following linear inequalities a > 1, b > 1, a > c, 
and b > c. The other condition is relaxed by minimizing d value. 


For the smooth variation of each matrix Af element due to V and a changes, the derivatives 
of each element with respect to‘ V and a are imposed by solving the following optimization 
problem. 

min W t [t ey eJ T (48) 

subject to 


1 

1 




1 

T 

9b 


1 


1 

T 

9 C 


01 


0 

T 

9d 


02 


0 

T T 

9 a ~ 9 c 


_03_ 


0 

< 3 - ^ 
1 


1 

0 


[1 <pi (p 2 (pz]gd < t, 


(49) 
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| ^(a^+^ini) | 

1 9V 1 

< e Vi , 

1 d(a^ 1 +(f> 1 n 1 ) I 
1 da 1 

^ £a± 

| 9(af 2 +0in 2 ) | 
1 9V 1 

< e V 2 , 

| 9(a\ 2 +(j>\n 2 ) , 

1 da. ' 


| 9(a% 1 +<f> 2 ni) | 
1 9V 1 

< e V 3 , 

I d(a 21 +<fi 2 ni) , 
I da 1 


1 9{a v 22 +<f> 2 n 2 ) | 
1 9V 1 

< e V 4 , 

| d{a 1 22 +(j> 2 n 2 ) , 
1 da 1 

^ ^ a4 

| dia^+foni) | 
1 9V 1 

< e v 5 , 

| d(a p 31 +(j> 3 n i) , 
1 da I 

< €a 5 

1 9(a^ 2 +rf> 3 n 2 ) | 
9V 

< e V 6 , 

1 9{al 2 +(j> 3 n 2 ) , 
da 

^ ^aQ 


( 50 ) 


where W t is a constant weighting vector to scale the parameter t, 6y, and e a . The derivatives 
are calculated using the first order finite differential method between grid points. 


After solving the optimization problem using the LMI solver “yalmip’ [11], the eigenvalues 
of the matrix Af at every grid point are shown in the right plot of Fig. 3. In Fig. 4 variation 
of (in of the matrix Af due to V and a is shown to verify the smoothness of the solution as 
compared to the particular solution calculated. 


For comparison of the nonlinear dynamics and the generated quasi-LPV model, the dy- 
namics are simulated for the initial condition [V a q 6} = [67 6 0 6], with control 
action changes such as the thrust doublet and the up-and-down elevator changes, respec- 
tively. The simulations results and the control command history are shown in Figs. 5 and 6. 
Note that the generated quasi-LPV model captures the time history of the original nonlinear 
dynamics. 


6 CONCLUSION 


In this paper, a new quasi-LPV representation method is proposed to reformulate a nonlinear 
mathematical model into a quasi-LPV form. The suggested method preserves the local and 
global stability of the original nonlinear system over the entire parameter space including 
the non-trim region. Nonlinear equations of motion for a general aircraft are formulated into 
matrix form which can be converted into quasi-LPV model using this method. For prac- 
tical use of quasi-LPV models, the nonlinear equations of motion are generally decoupled 
into the longitudinal and lateral-directional motions. The quasi-LPV representation of the 
longitudinal motion for the FASER aircraft has been developed, which has scheduling pa- 
rameters: velocity and angle of attack. The lateral-directional motion for the aircraft is also 
modeled as quasi-LPV representation with exogenous scheduling parameters V and a and 
state-scheduling parameters [3 and <fi. It has been shown that the longitudinal motion rep- 
resented by the quasi-LPV model preserves the stability of the original nonlinear dynamics 
at every grid point and represents the nonlinear dynamics over the entire parameter space. 
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APPENDIX A 


Matrix Form of Nonlinear Equation of Aircraft Motion 

After substitution of aerodynamic coefficients of Eqs. (36)-(41) into nonlinear force equations 
in Eqs.(33)-(35), the force equations are rewritten as 


Aforce\V Oi [3 p Q V B f orce [T S e &ra ^la &rf & if 6 r ~h foforce (51) 


where 


A 


force 


g g <l SS P C Yp qSbSpCy , 


0 0 


0 


2 mV 

C a Sf) -| qSc 


^{C x ,c a + c z ,s„) 


c„ 1 + - C Xq S a ) 

0 0 s a + ^Cy r -gjf (C x ,c a + c z ,s„) 


mV 


B 


force 


Bfi Bf 2 Bf 3 


qSbSpCYr 
2 mV 
S a Sf3 
Cp 

I qSbCp 

Ly ° l + 2 mV 2 L " Y r- 


(52) 

(53) 


B f i - 


CoiCp qSCp / | q \ 

m m ^ ^s e & Zs e a ) 

-J ir t Jk( C Z’.C«-Cx,.S«) 
_cys t _1 ^> (C x n+C z,S a ) 


mV 


mV 


(54) 


B f2 - 


qS_ 

mV 


CpV ( Cx Sra C a + Cz Sra S a ) + SpV C'Y Sr 
(C Z5r C a -C Xsra S a )/C 0 
. C 0 C Ysra - c a s 0 c x , - SpS a Cz. 


Cpv (C Xha C Q + Cz Sla S a ) + SgV Cy Sr 
(' Cz h C a - c XS[a s Q )/c 0 
C 0 C Ysia - c a s 0 c Xs ‘ - s 0 s a Cz Sia 


(55) 


Bf3 - 


and 


qS_ 

mV 


CpV(C Xtr C a + Cz Srf S a ) + SpVC Ytrf CpV(C Xsi C a + Cz tlf S a )+SpVC Ysif s 0 vc Y , 


°rf 


°rf 


If 




«lf 


(c ZSrf c a -c XSrf s a )/c 0 

CpC YSrf - C a S 0 C Xg f - Cz^ t SpS a 


°rf 


°rf 



\ ~ qSC e (C x C a + Cz S a ) 1 

ffl \ zv o or 1 Zj o or/ 


' CeSpSp - S e C a Cp + CeCpSaCp " 

h force 

Jk, - f (c z ,c a - c x ,s a ) 

+ 

y^(C a C e Cp + S e S a ) 


+ C Z ,S„)_ 


MCeSpCp + SpC e C a - C e C 0 S a Sp)_ 


(C Zsif C a -C Zsif S Q )/C 0 0 

c 0 c Ysif - C a S 0 C Xsif - Cz Slf s 0 s a C Ysr Cp 

(56) 


( 57 ) 
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Here, S a , C a , Sp, Cp, Sp, Cp, Sq, and Cq denote sin a, cos a, sin/3, and so on. 

The gravity term in Eq. (57) can be partitioned into 9 and 0 terms using polynominal 
expressions of cos and sin functions. In this report, 

cos(x ) =1 + b 1 x 2 + b 2 x 4 
sin(x) =a\x + a 2 x 3 + 03 a ; 5 

The gravity term can be written as: 


9vd9 

9v<j>9 


9 

0 


" S a Cp 

VC p 

vcp9 a< t > 


+ 

i r' i 

VCp a 

. y9ps 

ydsp . 


L T J 


~^SaSp_ 


where 


g V Q = — (cii + a 2 9 2 + a 3 9 4 )C a Cp + ( b\9 + b 2 9 3 )S a Cp (59) 

g v p = (01 + 02(1 + bi9 2 )(j) 2 + + a\b 2 9 4 )Sp + ((61 + & 2 0^)0 + b 2 (f> 3 )S a Cp (60) 

g a e = (ai + Qj 2 6~ + a^9 4 )S a + {b\9 + b 2 9 3 )C a (61) 

9a p = {h 0 + &20 3 + b\ 9 2 (f))C a (62) 

gpe = {he + b 2 9 3 )C a Sp - {he + b 2 9 3 )S a S/3 (63) 


gpp — (oi T 02(1 + b\9 2 )(j) 2 + ct30 4 + + aib 2 9 4 )C p — (&i0 + b 2 (f> 3 + b~l9 2 (j))S a Sp (64) 


The moment equation is written as: 


P 

Q 

r 


J 4 A m [ V a P p q r] T + J 1 B, m [T 5 e S ra ha hf hf h] T + J 1 h m (65) 


where 


A r 


0 

0 

qSbC h 

+ JzxQ — Jxyf 

JyzQ Jzz r 

Q 2V + JyyQ Jyz r 

0 

0 

0 

Jzz T J zxP 2V l q JyzP 

J ZX T* JxxP 

0 

0 

qSbC n/3 

q 2 V ^n.p + JxyP JyyQ 

J xxP JxyQ 

^§^rC nr + Jy Z p — J zx q 


u zx' 

( 66 ) 



b 

0 

0 


0 

0 

Cl, 

1 Ora 

Ci, 

5 la 

Cl , 

l6 rf 

Cl, 

s lf 

Ci Sr 

Brn = qS 

0 

c 

0 


0 

r< 

m s e 

r 

° m Sra 

r 

0 m ha 

r 

m S rf 

r 

^ m hf 

r 

m Sr 


_0 

0 

b_ 


0 

0 

c ns 

l °ra 

c ns 

US la 

r 

0 US rf 

r 

ns if 

r 


( 67 ) 
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and 


hm 


0 

qScC mo 

0 


( 68 ) 


Using Eqs (51) and (65), the full nonlinear equations of motion are written as: 


X* = A(V, a, 0, p, q, r, 9, <j))X st + B(V, a, (3)U + h(V, a, (5) 


Here, 


where 


and 


and 



I 3x3 03x3 03 x 2 


A r A 

force ■ rL g 

A(X st ) = 

03x3 J 1 03x2 


A m 03x2 


_02x3 02x3 hx2_ 


Atp 02x2_ 


0 0 0 0 C* -S# 


0 0 0 1 tan 6 tan 6 


A 


9vd9 9v<j>9 

VCp9 ae \ 

. v9bs vQpb . 



I 3x3 03x3 03x2 


B force 

B = 

03x3 J 1 03x2 


Bjn 


_02x3 02x3 hx2_ 


_ 02x7 _ 


I 3x3 O 3 X 3 O 3 X 2 


h force 

h = 

03x3 J 1 O 3 X 2 


h"m 


02x3 02x3 1 2x2 


02 x 1 


(69) 


(70) 


( 71 ) 


(72) 


(73) 


( 74 ) 
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APPENDIX B 


Matrix Form of Longitudinal Motion 

Consider longitudinal motion only. To extract the longitudinal motion from full equation of 
motion in Eq. (69), the state values related to lateral-directional motion are set to zero and 
the components ( J yx , J yz ) of inertia moment matrix J are set as zero because of symmetric 
configuration about the x-z plane. The longitudinal motion is rewritten as 


where 


V 


V 


T 

a 


a 

+ Bi on (V, a) 

Se 

Q 


g 


°a a u 

9 


9 - d 0 . 


fall _ 


+ ^hon(C, n) 


A 


Ion 


0 0 
0 0 
0 0 
0 0 


^(C Xq C a + C Zq S a ) 

1 + 2 ^(C Zq C a ~C Xq S a ) 

qSc 2 p 

2 VJyy 0m '3 

1 


g(— C a cos 6*o — S a sin 9 0 ) 
g(S a cos 9 0 — C a sin 9 0 ) /V 
0 
0 


(75) 


(76) 


Recall that the longitudinal control surfaces of this airplane are thrust, an elevator, ailerons 
and flaps. The matrix Bi and B 2 are for the conventional control surface (T, 5 e ) and for the 
ailerons and flaps, respectively. 

Bion = [B-i B 2 \ (77) 


where 


Bo = 



cy 

m 

f( C.'C'. + P'SJ 

Bi = 

s a 

mV 

5„) 

0 


qbc 


0 


0 

~{C X , 

m \ x s a 

C Q + 

c«.s a ) 

f (a,,c„ + 

qS p 

mV^ z Sa^ a 

c ZSa s a ) 

qS (p p 


qSc/^1 

j O 77 

J yy 

0 


qSc/^i 
t '- yl n 

0 


J yy 
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The nonlinear components of h are rewritten as 


h 


Ion 


~;(C Xo C a + C Zo S a ) + g(—C a sin# 0 + S a cosOo ) 
C Zo C a — C Xo S a ) + g(S a sind 0 + C a cos9q)/V 

J yy 


o 


Note that the aerodynamic stability derivatives are functions of angle of attack a. 

Matrix Form of Lateral-Directional Motion 

Under assumption that q is zero and V, a, and 9 are at a trim condition, lateral-directional 
motion can be extracted from Eq.(69). The lateral-directional motion can be rewritten as: 

Xiat — Aiat(V, a, /3, 9, 4>)X lat + Bi at {V , a, (3)Ui at (78) 


Here, 



i 

■ca 


Bp 

A-lat 

A 

x*- pr 

i Bi a t 

B pr 


_A ( j > _ 


_B(j)_ 


(79) 


where 


An = 


A 

-tipr 


Atj, — 

Bp = 


B pr 


Bfj, — 

and 


qSCpCy^ 


(Sc + ggKV,) (~Cc + WiCy,) f 


qSbCp 


mV 


J. 


-i 

sub 


qSbc h &£c, r *£ c lr o 

f ,2 


L qSbC np c np ^rC nr 

[0 1 tan 9 0] 

(^(CqY^-CcSqC^-SqScC^J 3(C^ - C a S„C x , 


qSbJ , 


-l 

sub 


Ua as 

c, 


l S a 


c L 


° a as 

C, 


S 0 S a C Z ) CpCy Sr 

Jas ' 


a 


n Sa 


a 




[0 0 0] 


hr 


a 


ns r 


(80) 


J sub 


J xx 

_ Jxz 


J X 

j zz 


(81) 
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APPENDIX C 


Characteristic Polynominal for Longitudinal Motion 

The characteristic polynominal of the matrix Af described in Eq. (46) is written as 

s 4 + as 3 + bs 2 + cs + d. 

The coefficients a, b , c and d are follows: 

a = g T a x 

where gj = [-(a^ + a 22 + 033 ) — ri\ — n 2 0] and x T = [1 (pi <p 2 03]. Recall that ri\ and 


n 2 are 

the elements of the nullspace A/". The coefficient b is 



b = g£x = [ 6 0 b\ b 2 b 3 \x 

where 

b 0 = 

( n P _l n P \ n P _ n P n P 4 . n P n P - n p a p - n p a p 
\ u 22 ' u UJ u, 33 u 32 u 23 ' u ll u 22 u 21 u 12 u 31 u 13 


61 = 

( a 22 A Q 33 ) n 1 — a^ 1 n 2 


b 2 = 

(033 + a 41 )n 2 - a 42 n, 


h = 

» V 

^23^2 ^13^1 ’ 


The coefficient c is 


c = g T c x = [c 0 ci c 2 c 3 \x 

where 

c 0 = 

( a 21 a i2 - «?i^ 2 ) + (<&<& - aV? 2 W 3 
+ («?1«23 - aP 2A) aP 32 ~ (<4 + ^S^Kl 


Cl = 

( a 32 a 23 ~ «33 a 2 2 )^l + («2l4s ~ a Il a 2 3 ) n 2 


c 2 = 

( a i 2 a 33 - a§ 2 ak)ni - (^a? 3 - a p u a p 3 )n 2 


c 3 = 

( a l3 a 22 - a 23 a l2 - < 4)^1 
+ ( a ll a 23 “ a 2 l a 13 - <* 24 )^ 2 - 

The coefficient d is 


d — gJjX — [d 0 d\ d 2 d^\x 

where 

do = 

( a ii a | 2 - a§ia? 2 )a §4 + (<&<& - a p 21 a p 2 )a p u 


d\ = 

p p P P 

^24^32^1 ^31^24^2 


d 2 = 

p p P P 

Q'31^'14^'2 ^32^ 14^1 


d 3 = 

(af 2 a? 4 - a? 2< 4)m + (a?^ - c4<> 2 . 
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Figure 1: The FASER 



Figure 2: The numbered equilibrium points and grid points 
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Figure 3: Eigenvalues of particular solution A p and the matrix A of the quasi-LPV model 
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Figure 4: Variation of an element of the matrix Af 
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Figure 5: Time responses due to the thrust level change 
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Figure 6: Time responses due to the elevator deflection change. 
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